2 



Nanoindentation of a circular sheet of bilayer graphene 

M. Neek-Amal^ and F. M. Peeters^ 
^ Department of Physics, Shahid Rajaei University, Lavizan, Tehran 16785-136, Iran 
Departement Fysica, Universiteit Antwerpen, Groenenborgerlaan 171, B-2020 Antwerpen, Belgium 

January 20, 2013 



Abstract 

Nanoindentation of bilayer graphene is studied using molecular dynamics simulations. We 
compared our simulation results with those from elasticity theory as based on the nonlinear Foppl- 
Hencky equations with rigid boundary condition. The force deflection values of bilayer graphene 
are compered to those of monolayer graphene. Young's modulus of bilayer graphene is estimated 
to be 0.8 TPa which is close to the value for graphite. Moreover, an almost flat bilayer membrane 
at low temperature under central load has a 14% smaller Young's modulus as compared to the 
one at room temperature. 

1 Introduction 

Graphene is an almost flat one-atom-thick layer of carbon atoms that are densely packed in a hon- 
eycomb crystal lattice. Most of previous studies concerned the electronic properties of graphene, 
but these two dimensional graphene membranes have also exceptional mechanical properties [1, 2, 3]. 
What makes graphene so exceptional is that it stays strong and stiff even up to a single atomic layer 
which differs from most materials whose mechanical properties substantially deteriorate as they are 
made thinner. These mechanical properties make graphene suitable for applications in e. g. pressure 
sensing. A nano-membrane made of graphene, the ultimate down-scaling limit of a nanoelectrome- 
chanical system, acts as a sensor that is predicted to detect ultra-small masses and being at the same 
time extremely robust against long-term wear. 

Recently, a nonlinear behavior of stress-strain response of monolayer graphene was studied by 
Lee et al [4] using atomic force microscopy (AFM). They measured Young's modulus and found 
that graphene is a strong material like diamond. Furthermore, linear force-displacement curves were 
measured in Refs. [5, 6, 7]. They obtained an effective spring constant for a micron size monolayer 
graphene sheet equal to 0.2 N/m [7] and for a suspended micron size sheet of multilayer graphene 
with thicknesses in the range 2-8 nm and found values in the range 1.0-5.0 N/m [5]. Moreover, 
monolayer graphene has a negative thermal expansion up to 900 K and its Young's modulus increases 
with temperature up to 900 K [8]. These anomalous properties are a consequence of the strong 
anharmonicity in graphene. Most of these properties can be explained by traditional elasticity theory 
and statistical physics [5]. 

The aim of the present paper is to investigate if two van der Waals coupled graphene layers, also 
called bilayer graphene, has similar improved mechanical properties as monolayer graphene. 

There are three common theoretical methods to study the indentation of a graphene sheet, fi- 
nite element method, continuum mechanics and molecular mechanics simulations [9, 10, 11]. These 
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methods lead to different results for the force deflection curve and for the graphene deformation pro- 
files [11]. The discrepancy originates from the predominant bond stretching mode predicted by the 
molecular mechanics model and a bending to stretching transition process under increasing deflection 
predicted by continuum mechanics. Nowadays various properties are measured on micron size circu- 
lar graphene [4] , while computational studies on nanoindentation of graphene typically are limited to 
nano size systems. In a typical experimental setup, the graphene sheet is indented by a tip from an 
atomic force microscope (AFM). 

An important issue in simulating nanoindentation of graphene is the type of interatomic potentials 
between the carbon atoms in the sheet and between the carbon atoms and the atoms of the tip. In 
recent studies on nanoindentation of graphene several models were used. Xu et al [11] used 1392 
atoms (a circular graphene sheet with radius 3.25 nm) in their atomic simulation and obtained the 
elasticity response of such a nano circular sheet of graphene under central load. The nanoindentation 
was realized by moving down the center of the sheet by hand. On the other hand, Hemmasizadeh 
et al [10] compared the results for the force deflection curve obtained from a continuum elasticity 
model [13] and finite element method calculations for two circular and hexagonal continuum plates 
of size 11 nm and showed that in the large deflection regime there are important differences between 
these two approaches. This difference roots in the effect of exerting a concentrated force so that the 
difference reduced to half when a spherical indenter with radius 10 A was used in the finite element 
method calculation. 

Medyanik et al [9] studied the nanoindentation of a hexagonal graphene sheet and multi-layered 
graphite sheets by employing a quasi-static formulation of the method of multi-scale boundary condi- 
tions. They showed that their method of multi-scaling gives good results in comparison to full domain 
atomistic simulations based on molecular mechanics simulations. They used repulsive interactions 
to model the spherical indenter. The number of atoms in the reduced domain was 2646 for single 
circular graphene. Depending on the size of the reduced domain, their results for multilayer graphene 
showed an increased force deflection dependence as compared to monolayer graphene. 

In this paper we study the mechanical properties of bilayer graphene for three different sizes 
and compare the mechanical response of bilayer graphene to the one of monolayer graphene. The 
Young's modulus of bilayer graphene is estimated by using the predictions from the theory of elasticity 
for a loaded plate in the large deflection regime. Furthermore, the temperature dependence of the 
nanoindentation of bilayer graphene is investigated. 

This paper is organized as follows. In Sec. 2 we will introduce the atomistic model, the simulation 
method and the the results from elasticity theory. Section 3 contains the numerical results and in 
Sec. 4 we will present our conclusion. 

2 Theory and model 

We have used classical atomistic molecular dynamics simulation (MD) to simulate the nanoindentation 
of a suspended sheet of bilayer graphene. The number of carbon atoms, i.e Nj,, varies from 35,688 to 
152,308 which are equivalent to surfaces with radius _R=12 nm and R=25 nm, respectively. A rigidly 
clamped boundary condition was imposed. In our simulations, the indenter consists of A^t=371 atoms 
in a fee structure (lattice constant equal to 3.92 A) and it is assumed rigid during our simulation. The 
shape of the indenter was chosen as a square-based pyramid. The area of the bottom surface (square) 
of the tip is 2.02 nm^ and the top one has an area of 10.24 nm?. Initially the coordinates of all atoms 
in each layer are put in a flat surface of a honey comb lattice with nearest neighbor distance equal 
to 0.142 nm while the upper layer is shifted along armchair direction by 0.142 nm above the bottom 



2 



layer. The separation between the two layers is 3.5 A. The initial velocities in each direction, were 
extracted from a Maxwell-Boltzman distribution for the given temperature. We simulated the system 
at room temperature 300 K and 20 K by employing a Nos'e-Hoover thermostat. The Brenner's bond- 
order potential [14, 15] was used for the carbon-carbon interaction and a Lennard-Jones potential 
(LJ) U{r) = 4£{{a/r)^^ — (cr/r)^} for the indenter-graphene interaction and the interaction between 
the two graphene layers. In the LJ potential, a is the distance at which the potential is zero and e is 
the depth of the potential well. The van der-Waals interaction between the two graphene layers was 
modeled by a LJ interaction with £c = 2.84 meV and ac =3.4 A. For the interaction between the tip 
and graphene, we used the LJ parameters for Pt atoms with ept=68.3 meV and apt =2.54 A [16]. For 
a two-component system, as studied here, the parameters for the mixed interaction between the two 
type of atoms can be estimated by the simple average (Jc-Pt = (fc + <ypt)/'^ and ec~pt = ^J£C ' £pt 
suggested by Steel et al [17]. To save computational time, we truncated the LJ potential at the cut-off 
distance of Tc = 3.5 a. Note that the LJ potential is a simple choice for modeling the interaction 
between two layers [18]. To obtain more accurate results one can use other potentials such as a Morse 
potential or other force fields [16]. At the start of our simulation, the position of the lowest atoms of 
the tip are located a few angstroms, i.e. Ri3.4 A above the upper graphene layer. After equilibrating 
the system during 50000 time steps, the indenter is pushed down slowly with (5=0.2 A in a time 
span of 5000At which is equivalent to a velocity of 8 m/s, where At=0.5 fs is the time step in our 
simulation. To avoid unphysical effects due to the time step, the indentation step -5- was chosen to 
be small with respect to the force cut-off length -rc- for the interatomic potential (e.g. to prevent a 
sudden change in the force, etc). 

The considered size of the system (~ nm) is larger than the deflection value (~ A), and in addition 
the thickness of the sheet is smaller than the amount of deflection. Therefore, nonlinear elasticity 
theory [19] for a circular flake in the large deflection limit along the z direction is applicable. The 
energetics of such a circular flake in the limit of large deflection is considered including both bending 
and stretching energies [19]. The condition of minimum energy for the flake yields the Foppl equation. 
The solution of the equation is obtained by using the Hencky transformation. The governing equations 
in planar-polar coordinates are as follows 

dr rdr ^ 2 dr ' 

dz F 

where r is the radial position, R is the radius of the circular plate as shown in Fig. 1, z{r) is the 
deflection at radial position r, t is the thickness of the plate, E is Young's modulus, cr,. is the radial 
stress of the flake and F is the concentrated load on the flake. We are interested in the situation of 
a rigidly clamped boundary condition or fixed boundary condition in the absence of residual stress, 
i.e., z = at r = 0. Expressions given by Eq. (1) are nonlinear equations, however they can be solved 
analytically in the special case of fixed boundary condition [20]. In general, the solution is given 
by [20] 

„ TvEt 1 , ,o 

where is a complicated function of the Poisson ratio, v. G{i') has almost a linear dependence 

on the Poisson's ratio of the desired system and varies from 0.9 for v ~ 0.0 to 0.65 for v ~ 0.5. 

C, is taken as the deflection of the bilayer graphene at r = 0, i.e. C, = z{Q). It is the difference 
between the center of mass of the central point of indented bilayer graphene and the first central 
non-indented equilibrium position of the bilayer (at r = 0)[19] (see Fig. 1). During indentation. 
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because the using of the LJ potential, the equihbrium distance between the lowest atoms of the tip 

and the central point of the bilaycr is of the order of cic-Pt- Surprisingly, our computer simulations 
confirm this behavior which we will discuss in the next section. The force-displacement curves have 
been measured recently by Lee et al [4] and they showed that it can be approximated by a simple 
polynomial function having a linear and a cubic term, 

F = aC + bC^. (3) 

When the bending stiffness is negligible and the load is small the force deflection can be approximated 
by the linear term while the second term dominates for large deflection. This behavior was studied 
both theoretically and experimentally in Refs. [4, 13, 21]. 

3 Numerical results 

Figure 1 gives a schematic view of the system showing the indenter and the clamped circular bilayer 
graphene and a snapshot of an atomistic indenter over a circular bilayer graphene with clamped 
boundary conditions. The z-component of the forces from the bilayer graphene atoms on the indenter 
are calculated by summing over the total reaction forces: 

Nt Nt 2 Nt Nt rr ^ 

F = T.T.F'^^ = 24eE{EE(-l)'2'=-l(-^)6'=i%}. (4) 

i=lj=l k=l 1=1 j=l 

Often in molecular dynamics simulations one approximates the above sums by including only the 
nearest neighbors in order to reduce the number of interactions which is accurate in the case of short 
range potentials. Regarding the cut-off distance (rc), only those bilayer atoms below the tip interact 
most strongly with the tip atoms, while outside this region the interaction strength decreases very 
fast. Therefore, in practice the sum over Ni, can be truncated and limited to the atoms below the 
tip. This is done by employing a neighbor list in our molecular dynamics simulation. The upper 
(bottom) layer atoms of the bilayer, below the tip region, interact with maximally Nt ~300 (250) tip 
atoms for small deflections and Nt ~320 (265) for large deflections. Moreover, the first term in Eq. 
(4), i.e. k=l in the parentheses, is the derivative of the attractive part of the LJ potential and the 
second term {k=2) is related to the derivative of the repulsive term. We will compare our obtained 
Young's modulus with those measured in experiments on graphene and graphite [4, 22] and compare 
our simulated force deflection curves of the clamped circular bilaycr graphene with those obtained for 
monolayer graphene [23]. Circular red data in Fig. 2 show the variation of the applied load at r = i? 
as a function of the deflection in the z-direction. The fluctuation in the data are larger than those 
in our previous study on monolayer graphene. The reason is that in bilayer graphene both layers 
vibrate (due to thermal fluctuations) almost independently. 

The plotted force was obtained as follows. After each 5000 time steps the tip is pushed down with 
0.2 A in order to induce the deflection During this time interval we let the system equilibrate 
and in the last 1000 time steps of these intervals we calculate F and obtain the mean value of F. 
Thus in a simulation with 10^ time steps we have 200 points in Fig. 2. The tip atoms are closer 
to the upper graphene layer and the other layer is, on the scale of the tip interaction potential, far 
from the tip. The latter layer is more free to vibrate and it has a smaller interaction with the tip 
atoms. Therefore, the contribution of the bottom layer to the summation of Eq. (4) is less than the 
contribution of the other layer. In our simulation wc did not observe any defect formation even for 
large deflections. Moreover, all the bilayer data of Figs. 2 and 5 show steps in the force-deflection 
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curves, that correspond to structural relaxation in the process of indentation. The reason is that the 
number of bilayer atoms which repel the tip are almost constant during the indentation steps. The 
reason that the bilayer is a discrete sheet and when it is indented by 2 A the number of tip neighbors 
(from bilayer atoms), which interact with the tip, are changed non-continuously. Such steps can be 
smoothed by time averaging. As we see the curve for R=12 nm is much smoother (than the two 
others in Fig. 2), since it was averaged over six simulations with different initial velocities. 

The numerical results in Fig. 2 arc fitted to the expression given by Eq. (3) and are shown by the 
solid curves on the circular red data. The values for the fitting parameters a and b are presented in 
Table 1. The results in Fig. 2 show clearly the size dependence of the force-displacement curves and 
its difference from similar results for monolayer graphene which are shown by triangular green data 
in each panel of Fig. 2. The behavior of the tip-graphene flake interaction even for a few angstroms 
displacement can be understood from our simulation given here. To obtain quantitative results, we 
start with Eqs. (2) and (3) to describe the dependence of parameter b on the radius for only large 
deflections. It is easy to obtain an analytical expression for the b value for a large circular sample [20] 
which is given by 

■^Et 1 

Note that the parameter -a- in Eq. (3) is equivalent to an effective spring constant for perpendicular 
indentation of the bilayer graphene when one assumes the bilayer graphene as an elastic membrane. 
Furthermore, as can be seen from Fig. 2, for small deflection the forces arc almost linear, we used a 
linear fit {Fi{() = a'C) for small deflections (C < 1 nm of the order of \faji>) and for large deflections 
(C > 1 nm) we used a cubic fit (i^c(C) = ^'C^ + C). The obtained values for the parameter a' and b' 
arc listed in Table 1. For R=12, 15 and 18 nm the parameter C is 11.3 , 8.0 and 6.5 nN, respectively, 
which are of the order of a. 

The fitted curves are also shown in Fig. 2 by dash-dotted curves. Figure 3 shows the obtained b 
and b' values as a function of the bilayer graphene circular size. The fitted parameter b' follows very 
well a,l/B? function and the agreement is now better than for the previous b-values. 

For monolayer graphene, we fitted the function F = amC + bmC^ to the results reported in our 
previous study [24]. The corresponding values for a„ and b^ are listed in Table 1. The fitted curves 
are shown by solid curves. Analogous as for bilayer graphene we also used separate fits for small and 
large deflection: Fi{Q = a'„^( and i^c(C) = ^'mC"^ + C", respectively, which are presented in Table 1. 
The corresponding fitted curves are shown by the dash-dotted curves in Fig. 2. For R=12, 15 and 18 
nm the parameter C" is 8.1 , 3.7 and 3.9. nN, respectively, which are of the order of a'. 

Furthermore, we fitted the values for b' to Eq. (5) and obtained TTEt/G{v) 4 x 353.5 N/m. 
This result is obtained from the function G(z^)=-0. 591^+0.94 obtained from Ref. [20] and using the 
value v ~ 0.25 which is a typical Poisson's ratio for monolayer graphene [21]. To obtain the Young 
modulus we need an appropriate value for the thickness of bilayer graphene. One can estimate t as 
the distance between the two layers plus the thickness of a monolayer. For monolayer graphene some 
experimentalists used the value 0.23 A by assuming Young's modulus of graphene to be the same 
as in graphite [2]. Others use different values for the thickness of monolayer graphene. Saitoh et al 
estimated 0.874 A which was obtained from the parameters of the Brenner potential [24]. This is an 
independent and more accurate value. Using this value the bilayer thickness becomes 4.5 Awhich leads 
to the Young modulus 0.8 TPa which is close to the values for graphite (i.e. 1.02 ib 0.03 TPa [22]) 
and those, found in recent experiments and found in theoretical calculations, for graphene (i.e. 1.0± 
0.1 TPa [4]). 

In order to investigate the effects of having two instead of a single graphene layer we calculated the 
ratio Jb^'^y'''- which we show in Fig. 4. For small indentation the fluctuations in our results are larger 

^monolayer 
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R(nm) 


a(N/m) 


6(10^^N/m^) 


a'(N/m) 


6'(10^^N/m^) 


12 


7.10±0.4 


20.20±0.5 


7.15±0.9 


23.3±0.6 


15 


3.80±0.5 


15.60±0.2 


7.17±1.3 


16.77±0.8 


18 


4.60 ±0.5 


10.90±0.3 


3.19±0.9 


11.99±1.1 


25 


5.50 ±0.3 


4.20 ±0.2 


5.30±0.9 


6.04±0.9 


R(nm) 


am(N/m) 




<(N/m) 




12 


5.21±0.08 


10.30±0.07 


3.41±0.1 


13.48±0.1 


15 


3.16±0.06 


7.07±0.01 


2.65±0.3 


10.07±0.3 


18 


3.06±0.04 


4.46±0.02 


1.86±0.1 


5.30±0.2 


25 


2.53±0.07 


2.70±0.05 


1.50±0.2 


4.01±0.2 


30 


1.30±0.09 


2.03±0.06 


0.65±0.1 


2.30±0.2 



Tabic 1: Values of a and b determined by fitting Eq. (3) to our simulated data given in Fig. 2. Values 
of a' and b' arc determined by fitting separately the linear and cubic equations with the small and 
large deflection regions, respectively. The parameters am and bm and also and b'^ for monolayer 
graphene which we determined earlier [23] are also shown. 



because both layers vibrate more freely, while in the large deflection regime the elastic membrane 
energy is dominant. Furthermore, because of the nature of the LJ potential, the repulsion of bilayer 
graphene with the tip is more or less two times larger than for single monolayer graphene (see Fig. 4) . 
As can be seen from Fig. 4, for small deflection where the Hooke's regime is dominant, we see a large 
difference between the spring constant (a parameters) of bilayer graphene and monolayer graphene. 
In the linear regime in Fig. 4 we have J^^^'^^v^ ^ where a},{= a) and am are the effective spring 

"monolayer O-m 

constant of bilayer and monolayer graphene. Figure 4 shows that in the small (^—region we have on the 
average Fuiayer < "^^Fmonolayer while in the large deflection regime Fuiayer - '^Fmonolayer ■ Note that 
the ratio ^ in the large deflection limit is not exactly 2, because the bottom layer of bilayer 

^monolayer 

graphene is far from the tip with respect to the upper layer so the number of interacting atoms in 
bilayer graphene to the tip at fixed deflection is not exactly twice those of monolayer graphene. Note 
that for small deflection (when the sheet is not stressed) shoulders (with small upwards amplitudes) 
are seen in both layers of the bilayer (particularly in the bottom layer because it is more free) around 
and below the tip which attract (they are close to the tip) the tip and cause negative forces even for 
small positive deflections (as measured from below the tip). This attraction dominates the repulsion 
term that is due to the other atoms below the tip. This unusual behavior is more clearly seen for 
larger systems. 

At low temperature (r=20 K) graphene has a smaller roughness and behaves as a fiat honeycomb 
lattice. One can compare the low temperature data in Fig. 5 (an almost flat bilayer graphene with 
i?=12 nm) to those for room temperature. At low temperature and for a fixed value of the displace- 
ment the forces are smaller than those at room temperature. We fitted Eq. (3) to the data and found 
ai = 2 Ah N/m and hi = 18.5 x W^N/rn^. Using Eq. (5) yields Young's modulus 0.69 TPa for low 
temperature which is 14% smaller than Young's modulus for bilayer graphene at room temperature. 
Such an increase of the Young's modulus with temperatures is unusual but it is similar to what has 
been reported recently by Zakharchenko et al [8] for monolayer graphene. This unusual behavior in 
both monolayer and bilayer graphene is the consequence of strong anharmonicity in graphene. 
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4 Conclusion 



In this study we showed that bilayer graphene hke monolayer graphene and other carbon nanostruc- 
tures has an exceptional stiffness. Young's modulus for bilayer graphene was calculated and found 
to be 0.8 TPa. The force-displacement result could be fitted to the function F = aC + b^^. We 
found that for given displacement the exerted force on bilayer graphene has to be about twice the 
one on monolayer graphene. The force-deflection result is found to be temperature dependent. At 
low temperature Young's modulus is found to be 14% smaller than at room temperature. 
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Figure 1: (Color online) Top panel: Schematic model of the suspended bilayer graphene with rigidly 
clamped boundary conditions. Bottom panel: A snapshot of an atomistic indenter over a circular 
bilayer graphene with R =6 nm and clamped boundary conditions. 
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-t^(nm) 

Figure 2: (Color online) Load force as a function of displacement for three different radii. The 
analytical expression Eq. (3) is fitted and shown by the solid curves. The dash-dotted curves 
refer to two separate fits for small and large deflections. For comparative purposes we show 
also the results of monolayer graphene and corresponding fits [23] . 
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Figure 3: (Color online) The parameters h and h' as a function of the circular size of bilayer graphene. 
The inset shows the parameters hm and h'^ as a function of the circular size of monolayer graphene. 
In both cases the solid curves show al/E? function fitted to the h' and h'^ data. For the definition 
of b, hm and 6', h'^ we refer to the text. (i.e. Eqs. (3) and (5), respectively) . 
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Figure 4: (Color online) The ratio between the forces in loaded bilayer graphene and loaded monolayer 
graphene for three values of the size of the membrane. 
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Figure 5: (Color online) Comparison between the force deflection response for rough (i.e. T=300K) 
and almost flat (i.e. T=20K) bilayer graphene for a membrane with radius R=12 nm. Solid curves 
are the fltted results using Eq. (3) with fltting parameters a = 7.1 N/m and b = 20.20 x 10^^ N/m^ 
for T=300 K and also ai = 2.45 N/m and k = 18.5 x lO^^iV/m^ for r=20 K. 
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